arXivrl 506.07136vl [cs.CV] 23 Jun 2015 


Segmentation of Three-dimensional Images with Parametric Active 

Surfaces and Topology Changes 

Heike Benninghoff* and Harald Garckd 


Abstract 

In this paper, we introduce a novel parametric method for 
segmentation of three-dimensional images. We consider a 
piecewise constant version of the Mumford-Shah and the 
Chan-Vese functionals and perform a region-based seg¬ 
mentation of 3D image data. An evolution law is derived 
from energy minimization problems which push the sur¬ 
faces to the boundaries of 3D objects in the image. We 
propose a parametric scheme which describes the evolu¬ 
tion of parametric surfaces. An efficient finite element 
scheme is proposed for a numerical approximation of the 
evolution equations. Since standard parametric meth¬ 
ods cannot handle topology changes automatically, an 
efficient method is presented to detect, identify and per¬ 
form changes in the topology of the surfaces. One main 
focus of this paper are the algorithmic details to handle 
topology changes like splitting and merging of surfaces 
and change of the genus of a surface. Different artificial 
images are studied to demonstrate the ability to detect 
the different types of topology changes. Finally, the para¬ 
metric method is applied to segmentation of medical 3D 
images. 

Keywords: Image segmentation, three-dimensional 
images, active surfaces, parametric method, topology 
changes, Mumford-Shah, Chan-Vese, finite element ap¬ 
proximation 


1 Introduction 


One major challenge in image processing is the au¬ 
tonomous detection of objects in images and the seg¬ 
mentation of the objects from each other and from their 
environment. 

A very popular approach for image segmentation is the 
active contour method 20 , 17 . In the case of classical 


two-dimensional images, one or more curves, called con¬ 
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tours, evolve in the two-dimensional image domain and 
stop locally at edges or region boundaries. The motion is 
described by evolution equations which aim to minimize a 
certain energy functional. The energies typically contain 
length terms to control the smoothness of the contours 
(internal energies) and terms which push the contours to 
the desired region boundaries or to edges in the image 
(external energies). 

Two main classes of approaches can be distinguished: 
The first class are edge-based methods where regions are 
identified by their boundaries where the image intensity 


function rapidly changes 20 , 23 , 13 . The second class 


are region-based methods where the regions are charac¬ 
terized by the mean gray value or mean color, or by the 


texture or some other grouping 28 , 32 , 16 , 38 


Region-based active contours methods can also be ap¬ 
plied on images with so-called weak edges, i.e. edges with 


only small changes in the image intensity function 16 


and on images which contain regions which are groups 
of smaller objects [3]. Furthermore, the method is less 
sensitive to noise. If images with high noise have to be 
segmented, gradient-based approaches may get trapped 
at locations where the noise is high and the contours may 
not detect the real objects in the image. 

In this paper, we study volumetric, i.e. three-dimen¬ 
sional, images given by a scalar or vector-valued image 
function uq : D -A where D C M 3 is an open and 

bounded image domain. Real images are often defined 
on a set of N x x Ny x N z voxels (=volume pixels), where 
uo is locally constant on each voxel. 3D images may be 
reconstructed by certain 3D imaging procedures like com¬ 
puted tomography (CT) or magnetic resonance imaging 
(MRI), cf. 39 , pi. 


3D image segmentation aims at dividing a given im¬ 
age in connected regions, representing 3D objects in the 
image or their environment in the image domain 12. For 
3D image segmentation, the boundaries of the regions or 
objects have to be detected. These boundaries can be 
represented by a set of two-dimensional surfaces. 

The active contours concept 20 can be extended to 
the three-dimensional case. For 3D images, we can con¬ 
sider time-dependent two-dimensional surfaces (active 
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surfaces ) and evolution laws for the surfaces which at¬ 
tract them to region boundaries. In particular, we will 
study extensions of the Mumford-Shah model 28 and 
the Chan-Vese model [16 to the three-dimensional case. 

belong to the first works, where 
is extended to volumet- 


The articles IT , 18 


20 


the active contours model 
ric image data. There, an analytical framework is intro¬ 
duced for 3D deformable surfaces. For practical compu¬ 
tations however, the authors suggest to replace a given 
3D image by a sequence of 2D images and to apply the 
2D active contour model on each single image, followed 
by a 3D reconstruction of the surface. 

The geodesic active contours model 13 is a popular 
edge-based method. An extension of the geodesic active 
contours model to 3D image segmentation is proposed 
. The level set method 


14 


29 


is used to describe the 


40 , where 2D 


surface implicitly. 

The level set method is also used in 
and 3D active contour models are presented and applied 
on medical images. However, practical results are only 
shown for 2D images. A detailed literature study on 
3D brain cortex segmentation is given in 


22 


on segmentation of medical X-ray computed tomography 
and magnetic resonance images is presented in 35]. 

In 27], a combination of edge-based and region-based 
segmentation methods is proposed. Both explicit (tri¬ 
angulated surfaces) and implicit (level set) methods are 
implemented. For explicit methods a constant global 
topology is assumed. For images which require topology 
changes, only the level set method is applied. 

The Chan-Vese model 16 is used for 2D and 3D medi¬ 


cal applications in 33 to perform heart segmentation us¬ 


ing an iterative version of the Chan-Vese algorithm. The 
level set method with a finite difference scheme is used 
to solve the segmentation problem numerically. The level 
set method is also applied in |1 and [36] for active sur¬ 
faces. Applications using 3D medical data (i.a. lung and 
heart segmentation) are considered. In 24], the level-set 
method is applied for 3D cell membrane segmentation. 
An approach for tracking cells in 4D images (3D data + 


time) has been developed in 25 


Also some parametric approaches for surface evolution 
exist in literature: In 11 , a program called ”The Surface 
Evolver” is presented which computes evolving triangu¬ 
lated surfaces, where the evolution is driven by energy 
minimization problems with possible constraints. The 
program is able to perform topology changes like split¬ 
ting if this is instructed by the user. 

Another explicit method for evolving surfaces by using 
triangulated surfaces is proposed in 


12 


. Techniques are 
introduced for mesh quality improvement and for topol¬ 
ogy changes. Splitting is done by a combination of remov¬ 


ing degenerate elements and a certain mesh separation 
method based on duplication and separation of nodes. 
Surfaces are thus split if they become locally too thin. 
Merging is detected by searching for edges which are too 
close. 

Finite element approaches for surface evolutions are 
pursued in |'6' . There, the authors do not consider image 
segmentation applications, but surface diffusion. Several 
mesh quality routines like mesh regularization (keeping 
all angles of simplices at a node of the same size), time 
step control, refine/coarsening routines and angle width 
control are proposed. 

In this paper, we present a novel parametric approach 
for 3D image segmentation. We present a scheme for im¬ 
age segmentation describing the evolution of parametric 
surfaces. For the numerical approximation, the smooth 
surfaces are replaced by triangulated surfaces. We make 
use of a parametric finite element scheme based on 8|. 
There, a scheme is proposed for surface diffusion, (in¬ 
verse) mean curvature flow and non-linear flows. We use 
and apply this scheme to image segmentation with mul¬ 
tiple phases and regions. 

Our method also allows for topology changes which 
have not been addressed in [8]. We efficiently detect 
topology changes and perform modifications of the sur¬ 
face triangulations. In |l0], we considered segmenta¬ 
tion of two-dimensional images, and used and extended a 
method to detect topology changes 26 to handle a vari¬ 


ety of topology changes of curves. Topology changes in¬ 
volving surfaces are more complex compared to topology 
changes involving curves. For example, if a curve splits 
up in two subcurves, the discretization has to be modi¬ 
fied only at two points, see 10 for details. If a surface is 


split up in two subsurfaces, many triangles are located in 
a small volume. The pure detection of such a splitting is 
quite simple; the idea of an auxiliary background grid 26 


as used for curves can be extended to topology changes of 
surfaces. However, the modifications of the surface trian¬ 
gulation are not as straight-forward as for curves. In case 
of splitting, we will delete the involved triangles near the 
splitting point resulting in two surfaces with intermedi¬ 
ate holes. Then, we will close the intermediate holes by 
creating new triangles. Apart from splitting and merg¬ 
ing, further topology changes can occur for surfaces: An 
increase or decrease of the genus of a surface can also oc¬ 
cur, for example when a sphere evolves to a torus or vice 
versa. In summary, the execution of topology changes is 
the main additional challenge of 3D image segmentation 
with parametric surfaces. Therefore, one main focus of 
this paper is the detection, identification and execution 
of topology changes. 

The remaining part of this paper is structured as fol- 


2 
























lows. In Section |2j we present a region-based active sur¬ 
face model where we extend the Mumford-Shah model 
and the Chan-Vese model 


28 


16 to 3D images. We 


present an efficient parametric scheme for the evolving 
surfaces. Also multiple phases can be handled. The main 
part of this paper is the numerical approximation of our 
scheme and the handling of topology changes which is de¬ 
scribed in Section[3] A finite element scheme is presented, 
a corresponding linear equation is derived and some com¬ 
putational details are given including mesh quality as¬ 
pects and time step control. The detection, identifica¬ 
tion and execution of topology changes is described in 
detail including a description how to modify the triangu¬ 
lations after a topology change has been detected. In Sec¬ 
tion [4| we present results from segmentation of artificial 
test images and real medical images. We demonstrate 
the different topology changes which can occur during 
the evolution of surfaces. A final conclusion is drawn in 
Section [5j 


The first term in 0 penalizes the area of the surfaces, 
the second term does not allow u to change much in fi\T, 
and the third term requests that u is a good approxima¬ 
tion of Uq. 

We first consider two-phase image segmentation, we 
consider one closed, orientable surface T separating two 
disjoint regions Di and U 2 such that Q = fli UTU^. 
We assume that T is oriented by a unit normal vector 
field V pointing from U 2 to ^i- 

Furthermore, we consider a piecewise constant version 
of the Mumford-Shah functional: We search for a surface 
T and for an approximation u : U -> M of Uq which is 
piecewise constant in each region, i.e. u\Q k — c^^k — 1,2, 
such that 


£(r,ci,c 2 ) = 

= cr\r\+\( f (uq — C 1) 2 da; + [ (it 0 - c 2 ) 2 drc ] (2) 
rq J rq J 


2 Segmentation of three- 
dimensional images 

2.1 Region-based Active Surfaces 

We perform image segmentation by active surfaces, the 
surface-analogue to active contours 201, 17 . The idea of 
active surfaces is to let surfaces T(t), t G [0,T], evolve in 
time such that a certain energy functional is minimized. 
From the minimization problem, one can derive an evo¬ 
lution law such that the surfaces of T(t) are attracted to 
the region boundaries in the given image. 

Let Q C M 3 be open and bounded. We first consider a 
scalar image function uo : £2 —)> M. 

In this paper, we restrict on region-based methods 
for segmentation of 3D images because of advantages 
of region-based approaches compared to edge-based ap¬ 
proaches (cf. Section [I]). 

The Mumford-Shah [28 method for 3D images aims 
at finding a set of two-dimensional surfaces T = Ti U 
... U Tjsf c and a piecewise smooth function u : ft M 
with possible discontinuities across T approximating the 
original image uq. The energy to be minimized is 


E ms (u, r) = a|r| + [ || Vu|| 2 dx + X f {uo- U ) 2 da:, 

Jn \r Jn 

(i) 

where cr, A > 0 are weighting parameters and |F| denotes 
the total area of the surfaces belonging to T. (For non¬ 
smooth surfaces, we identify \T\ with the two-dimensional 
Hausdorff measure of T CM 3 .) 


is minimized. 

Similarly, the Chan-Vese functional 16 can be ex¬ 
tended to 3D images. The energy to be minimized is 

£(r,ci,c 2 ) «cr|r| +/x / ldx + \i / ( u 0 -ci) 2 dx 
J rq J rq 


+ As 


/ (uo - C 2 ) 2 
J rq 


dx , 


( 3 ) 


where cr, Ai, A 2 > 0, /i > 0 are weighting parameters. For 
fi = 0 and Ai = A 2 = A, this is the functional §. 

The energy defined in 0 depends on the surface T and 
on the image approximation u given by the coefficients c \, 
C 2 . For minimizing 0, we perform a two-step approach: 

First, we fix the surface T and consider variations in 
the coefficients <q, C 2 . Using the theory of calculus of 
variations we obtain the mean of the image function in 
fife for Cfe, k = 1,2: 


Cfe 


fn k da: 

fn k 1 d ' c ‘ 


( 4 ) 


Then, we fix c\ and C 2 and consider small variations of the 
surface T by smooth surfaces T(t) C U, t G (—e, e), with 
r(0) = T. Let Qi(t) and fU(£) be the regions separated 
by r(£). We define 


f(x,c 1 ,c 2 



{uo{x) - Cl) 2 , 
(u 0 (x) - c 2 ) 2 , 


if x G Qi (t), 
if x G u 2 (t), 


( 5 ) 


which is defined for a.e. x G U. 
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By using a transport theorem, we obtain 
d 


d t 


E(T(t),c 1 ,c 2 ) = 


t =0 

d 
d t 


[a ldA + \ f(x,c 1 ,c 2 ,t) 

t =o \ Jr{t) Jn 


dx 


= —a I k V n d A+ 


L 

A J ((m 0 - Cl ) 2 - (m 0 - c 2 ) 2 ) V n d A 

- J {(tk + F)' 


Vn dA 


where d A is the area element, V n is the normal velocity, 
k the mean curvature and F is an external force given by 

F(x ) = A ((u 0 (x) - Cl ) 2 - («o(x) - c 2 ) 2 ) ,fer. (6) 


surfaces Ti(t), t E [0,T], i = 1,..., TV#, which separate 
three-dimensional regions Qk(t)i k = 1,..., JVr. We as¬ 
sume that the surfaces are compact and oriented by unit 
normal vector fields z7^(., t) pointing from to 

£4+(i) (t), where fc ± (i) e {1,.. .,N R }. 

Let Xi(.,t) : Ti —)> M 3 , i = 1,..., TV's, be a smooth 
parameterization of r^(t), where is a two-dimensional 
reference manifold, for example the sphere = S 2 C 
M 3 . The normal velocity of I\(£) can be expressed as 
(Vn)i — ( %i)t - Vi- 

An approximation of the image intensity function uq 
is given by the piecewise constant function u(.,t) = 
Ek=iCk(t)xn k (t), where Xn k (t) is the characteristic func¬ 
tion of and c k (t) is the mean of uq in Qk(t)- 

For each surface, we define the external forcing term 

Fi(.,t) = A ((u 0 - c ft+(i) 0)) 2 - (w 0 - c fe -(i)0)) 2 ) (9) 


The fastest decrease of the energy is obtained for 


V n » ok + F. 


( 7 ) 


Also multichannel images with a vector-valued image 
function —>> can be handled. This involves 

vector-valued coefficients Cfc, k = 1, 2, and a modification 
of the external force to, for example, 


d 

F(x) - A,; [((ito)i(f) - (ci)i) 2 - ((uo)i(x) - (c 2 )i) 2 

2=1 


(8) 

where the subscript i denotes the i-th component of a 
vector, i = 1,..., d. For computation of the coefficients, 
each component of c k is set to the mean of the corre¬ 
sponding component of Uq in the region k = 1,2. 

In principle, also spaces like the HSV (hue, satura¬ 
tion, value) or CB (chromaticity, brightness) space can 
be used [ 2 ], [l5], [371. In these cases, the image function 
has values on certain submanifolds of M d . In 10 , we 
proposed a method to segment 2D images using the color 
space HSV and CB. The method can be transferred also 
to the 3D case. In many practical applications however, 
for example medical 3D image data generated by com¬ 
puted tomography (CT) or magnetic resonance imaging 
(MRT), the image function is often scalar-valued (cf. for 
example the lung image database of The Cancer Imaging 
Archive (TCIA) |3l], [|, [30 ). 


2.2 Parametric and Multiphase Formu¬ 
lation 

Equation 0 can be rewritten using a parametric ap¬ 
proach to describe the time-dependent surfaces. Fur¬ 
ther, we now consider a more general setup of multiple 


and obtain the following scheme for the surfaces: Find 
Xi( ., t) : Ti —M 3 and /^(., t) : —>> M, i m 1,..., Ns, 

satisfying 

(xi) t - Vi = crKi + Fi, (10a) 

Ap Xi = KiVi. (10b) 


Equation (10a) is a parametric formulation of (|7|) for mul¬ 
tiple regions. Equation (10b) relates the parametrization 
Xi and the curvature Ki, see e.g. 19 . The symbol Ap 


denotes the Laplace-Beltrami operator. Here, we use a 
small abuse of notation, i.e. we consider Ki and as 
functions defined on Tj, i.e. we identify Ki with o Xi 
and Vi with Vi o xi, i = 1,..., Ns- 


3 Numerical approximation 


3.1 Finite Element Approximation 

We introduce a finite element approximation for the 
scheme ( [lQ] ) which is based on a scheme developed in [8], 
where geometric flows of two-dimensional surfaces are 
considered. We extend the ideas to solve schemes like 


(10) which arise in image segmentation applications. 

Let 0 = to < t\ < ... < tM = T be a decomposition of 
the time interval into possibly variable time steps r m = 
t m + 1 — t m for m = 0,..., M — 1. 

Let Ns denote the number of surfaces and Nr denote 
the number of regions. Let the smooth surface r^(t m ), 
i = 1,..., Ns, be approximated by a polyhedral surface 
r™ of the form 

Ni,F 

r r = U (u) 

3 = 1 
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where crjj, j = 1,..., TV^, are disjoint, open simplices 
(also called faces) with vertices gjj, j = 1,..., TV^y. Fur¬ 
ther, let h := max i= i v .. ? iv s j=i,...,Ar i)i rdiani((j^) be the 
maximum diameter of a simplex of the triangulated sur¬ 
faces. The diameter diam(ojj) is defined as the maxi¬ 
mum distance between two points of a™-. 

Let X™ be a parameterization of T™ and let fi™, 
k = 1,...,7V#, denote the open, disjoint subsets of Yl 
separated by T™, i = 1,..., TV's. Thus, Yl™ is an approx¬ 
imation of Clk(tm) for k = 1,..., Nr. 

The new surfaces r^ +1 are parameterized over Y™. 
Therefore, we define the following finite element spaces 

W{Y m ) ^{{r ]l ,...,r lNs ) eC(TT,m x...x 
C(T% S ,R) is linear, 

Vi = 1,, iYs, j = 1,..., V,f} , (12a) 

Z(r m ) :={(7/ 1 ,...,rf JVs )eC(r-,K 3 ) x ...x 
C(r)O s ,R 3 ) :fh\*~ is linear, 

Vi = 1,..., A/g, j = 1,..., N i}F } . (12b) 


where Qij t , l = 1,2,3, are the vertices of crjj, \a™j I = 

§ll(Cj2 “ Cil) x (Cfs “ Cn)ll is the area of a Z and 

We assume that the vertices {/])!), , j = 1,..., 

7V^ ?j p, are ordered such that the unit normal V 7 ? 1 at T ™ is 
given by 



:= iT 171 - 

hj 


(gy - CtJ x (Cia - CtJ 

ll(Ci 2 -C;J x (Ci3-Ch)ll 


(16) 


points from to fiJJ* +(i) . 

We propose the following finite element scheme ap¬ 
proximating the scheme (10): Let T° be a union of poly¬ 
hedral surfaces approximating T(0) and let X° G V(r°) 
be the identity function on T°. Find X m+1 G ]3(r m ) and 
ft m+1 G VF(T m ), m = 0,1,..., M — 1, such that 


1 _ 


-, X F™)^- ( T( K m + 1 ,x)J 




h 

m 5 


v x e w(P"), 


(17a) 


The spaces W(Y rn ) and ]/(r m ) thus consist of scalar or 
vector-valued, piecewise linear functions defined on Y m . 

A basis of VF(r m ) is given by functions xTj := 
((xS)i. • • - 0C>s) e lF(r m ), where 

(x™M^i)=5ikS j i (13) 

for i, fc = 1,..., N s , j = 1,..., A*,y, l = 1,..., N k y. 

Note, that depending whether the domain of definition 
is r m_1 or T m , we can interpret X 171 as a different func¬ 
tion. In particular we have for m > 1, X 171 G y(r m_1 ), 
and for m > 0, G ]3(r m ) is the identity defined on 

pm 

For scalar functions u = (ui ,..., un s ), v = (t>i, • • •, 
vn s ) G L 2 (r5 n ,M) x ... x L 2 (r^,M) and for vector¬ 
valued functions u = (txi, ..., un s ), v = (t’l,... , Vjv s ) G 
L 2 (r5 n ,M 3 ) x ... x L 2 (r^ s ,M 3 ), we introduce the L 2 - 
inner product over the current polyhedral surface Y 171 as 
follows: 


(u,v 



u . v d A 



Ui. Vi d A. 


(14) 


If u, v are piecewise continuous with possible jumps 
across the edges of ojj, i = 1,..., TVs, j = 1,..., W,f, 
the mass lumped inner product is defined as 

N s Ni, F 3 

(«w)m ; = o E E ( 15 ) 

2=1 j = 1 1 = 1 


( K m+1 i^,ff) h m + (V s X m+ \V s ff) m = 0, 

VrjeV(r m ). (17b) 

Here, F m = ( F {",..., Fff s ) is defined by 

prm) ;=A ((Mxrm)) - ^ +(i) ) 2 + 

~(u 0 (xr( :«&•))-#- w ) 2 )> 

for i = 1 ,... , TV's- and j = 1 ,... , TV^y and c™ is set to 
the mean of uo in Yl™ for k = 1,..., Nr. 

We further introduce a weighted normal defined at the 
nodes = Q^j ^Y™ by setting 


-m ( -*rn \ -*m __ 

^i 


1 


E 




(18) 






'T'm ._ 


where for i = l,...,A(g and j = 1, 

{<i ■ Ci e <[} and A-. := 

Further, we set ^(gfl-) := ^ := 1| and 

= «, • • • ,C3% g ) and> = (<,..., 

As in we make a very mild assumption on the 

triangulations: 

(A) For m = 0,..., M, we assume that |crjj| > 0 for 
all i — 1,..., Ns and j = 1,..., A^ ?j p and for m = 
0,..., M — 1, we assume that dim span 
3. 
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The assumption (A) is only violated in very rare cases. 
For closed surfaces without self intersections, it always 
holds. 

Theorem 1. Let the assumption ( A ) hold. Then there 
exists a unique solution G E(r m ) x 

W( T m ) to the system ( pT[ ). 

Proof. (See also [ 8 ].) Since the system is linear, it is 
sufficient to show uniqueness. Therefore, we consider the 
following scheme: Find X G H(r m ) and k G W(T m ) such 
that 


- — (X,x0m + ^(«,X)^ = 0, Vx 6 ^(r), (19a) 

'T~m 

(Ki7 m ,ff)t + (V s x,v s ff) m = o, v?? e Z(r m ), (19b) 

holds. Testing < |19a[ ) with X = k and ( |19b| ) with fj = X 
leads to 


aT m (K, + (VgX, V s X) m = 0. 


( 20 ) 


It follows that Kij = 0 and Xij = Ci G R 3 for i = 
1 ,..., Ns , j = 1,..., Niy. Inserting n = 0 and X = 
C = (Ci,..., Cat s ) in (19a) results in 

(C, X i7 n ) h m = 0, Vx e iT(r m ). (21) 


Choosing x = xY!j (th e standard basis) and using (18), 
the definition of cJ-j, results in 

=0, Vi = 1,... ,7Vs, j = 1 ,..., (22) 


Finally, using the assumption (^4), it follows that Ci — 0 
for each i = 1 ,..., Ns- □ □ 


3.2 Solution of the Discrete System 

We define SX m+1 := X m+1 -X m . As SX m+1 and P m+1 
are uniquely given by their values at the nodes gjj, we 
consider them as elements in (M 3 )^ and R N , respectively, 
where N = Y^i=i Niy- We introduce the matrices M m G 
R NxN , 7V m G ( R 3 ) NxN and A m G ( r.3xs^nxn by 

/ ^ ••• o 

Afrn * = ; *. ; 5 

V 0 ... M% s 

Pa ••• 0 \ 

N m := : : , 

V 0 ... yy / 


: — 


A 


0 


A Ns 


where M' m e rW,vxjv 1iV} y g (rS^.vxjv,^ ^ g 
(R 3x3 ) iV, ’ v ' xA ? .,v s = 1, ■ ■ ■, Ns- Their entries are defined 
by 


[Mt, 

\K 

X 


\kl 


\kl 


\kl 


= (x™k,x™ir m , 

= <x&,xmO£,, 

= (vy k ,vy)mid3, 


(23) 


with i = 1 ,..., Ns , fc, Z = 1,..., TVyy. Here, M 3 denotes 
the identity matrix in M 3x3 . Further, we introduce 6 m = 
• • • ; fc m s ) e RjV defined by 

[bL]k := (FT,X a) 


m \h 
m') 


— 1 , . . . , TVg , k = 1 , . . . , TVj y. 

(24) 

The scheme (pT| ) can be rewritten to the following prob¬ 
lem: Let T cr be a polyhedral approximation of T(0) 
and let X° = (X°,...,X° Ns ) G (R 3 ) N with X 3 = 
(X? A ,...,xy v ) such that xr are the coordinates of 
the vertices of T^ for i = 1,..., Ns , j — 1,..., iV^y. For 
m = 0,..., M - 1 find SX m+1 G (M 3 )^ and ^ m+1 G R N 
such that 


(JT m 


—N t 

m 

Am, 


^m+l 

(LA m+1 


Nn^rn 

-A m X m 


(25) 

Applying a Schur complement approach, we can trans¬ 
form this system to 


„ra+l _ 


1 


= -M! 


-1 


—yr<5x m+i - k 


(26a 


- NrnMrr} Nm + A„ 


^m +1 = 


1 

<7 


+ -NmM^bm- 


(26b) 


Since the system matrix in (26b) is symmetric and pos 


itive definite under the assumption (A), there exists a 
unique solution. 


The linear system (26b) can be solved with an itera¬ 


tive solver, for example, with the method of conjugate 
gradients with possible preconditioning, or with a di¬ 
rect solver for sparse matrices. For the experiments and 
examples presented in Section [4] of this paper, we use 
a MATLAB built-in routine, a direct solver for sparse 
systems. Even for two-dimensional problems, which re¬ 
sult from the evolution of two-dimensional surfaces, the 
sparse direct solver is very efficient from a computational 
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3.3 Topology Changes 


Parametric methods cannot handle topology changes au¬ 
tomatically in contrast to other numerical methods like 
the level set method. During the evolution of surfaces 
singularities can occur like a pinch-off, see ©> i . In or¬ 
der to proceed after a pinch-off, the surface has to be split 
in two single surfaces. Other possible topology changes, 
that we will consider here, are merging of two surfaces 
and change of the genus (occurs for example during an 
evolution of a torus to a sphere or vice versa). 


In 26 , an algorithm is proposed to efficiently detect 
splitting and merging of evolving curves in M 2 , see also 


[5]. In 10 , we used and extended this algorithm to detect 
topology changes in 2D images. In this paper, we want to 
adapt this approach to the 3D case. We aim at detecting 
topology changes which could occur during the evolution 
of surfaces. Having found the location where a topology 
change occurs, we propose a method how to modify the 
triangulations. 


3.3.1 Detection of a Topology Change 


In 26 , a virtual, auxiliary 2D background grid is con¬ 


structed which covers a two-dimensional domain, and 
topology changes of curves are detected if node points 
from different curves or different parts of one curve are 
located in one array of the background grid. 

Motivated by this method for evolving curves, we pro¬ 
pose the following method to detect topology changes of 
evolving surfaces. The basic idea is the use of a uniform 
3D grid of cubes. A topology change may occur, if a large 
number of nodes or if nodes of different surfaces or dif¬ 
ferent parts of one surface (with opposite normal vector) 
are located in one cube. 

In detail, to detect a change in topology, we construct 
a uniform 3D background grid which covers the image 
domain D. In the following, we assume that D is a cuboid. 
If the 3D image uo is not given on a cuboid volume, we 
consider a cuboid which contains D. 

Let D [^min^max] X [|/min ? 1/max] X [^minj ^max] ^ M 

be the image domain containing in particular the surfaces 
r™, i = 1, ..., Ns- We consider a grid dividing Q in 
a set of many small cubes of edge width a E M. Let 
the grid consist of N x x N y x N z cubes, where N x = 

Ceil((x m ax %min)/N), Ny Ceil((^ max 2/min)/^) £md 

N z = Ceil((z max ^min)/^)* 

We now perform one loop over all surfaces and nodes 
X™j, i = 1 ,..., N s , j = 1,..., Niy. If a node is the 
first node which is detected to lie in a certain cube, we 
create a new list for that cube, where we store the index 
pair (z, j). If another node has already been identified to 


be located inside that cube, we add the index pair (z,j) 
to the existing list. 

If there is a large number of nodes located in a cube, 
i.e. more than Ahetect nodes, a topology change likely 
occurs, and the cube is stored in a list for possible topol¬ 
ogy changes. It is also possible that the node density is 
only locally very high at this location, but no topology 
change happens. 

If there are less than detect nodes in the current list, 
we compare the surface index and the direction of the 
weighted normal vector of the current node with those of 
the nodes already stored in the list. If two different sur¬ 
face indices i\ and 12 occur or if two nodes with (nearly) 
opposite weighted normal vector are located in one cube, 
a topology change likely happens. The corresponding 
cube is accordingly stored in a list for topology changes. 

After having considered all nodes, the cubes marked 
for topology changes are considered one by one. If a 
topology change is identified, the surface triangulation is 
accordingly changed. Thus, by successively considering 
all marked cubes, more than one topology change can 
be executed in one time step. The list can be optionally 
sorted such that the cube with the largest number of 
nodes is considered first. 

The detection of topology changes is very efficient from 
a computational view. The effort is of order 0(N ), where 
N is the total number of node points. For comparison, 
a simple approach, where all possible pairs of two nodes 
are considered and where the distance between two nodes 
is computed to detect a topology change, would result in 
a computational effort of order 0(N 2 ). 

Similar as described in |10], the grid size a can be 
adaptively set, for example dependent on the speed of 
the evolving surfaces. Therefore the method to detect 
topology changes is both efficient and robust. 


3.3.2 Identification of the Topology Change 


For identifying which kind of topology change occurs, the 
nodes of the affected cube, i.e. the current cube of the 
sorted list, and the nodes of up to 26 neighbor cubes (in 
total up to 27 cubes, i.e. 3x3x3 cubes) are considered. 
Let S = {j 1 ,..., j nc } denote the index set of the nodes 
and let Xj, j E S, denote the coordinates of the nodes 
located in the cubes. Further, let ujj, j E S , denote the 


corresponding weighted normal vectors at Xj , recall their 
definition in (18). For the ease of illustration, we omit 
the time dependency (time index m) in the notation. 

The different topology changes are distinguished by 
considering the weighted normal vectors uij. The idea is 
that in case of merging and in case of increasing genus, 
their are two main group of nodes which can be found 
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by considering their normal vectors. For splitting and 
decrease of genus, there are more than two main direc¬ 
tions. 

The node with index j = ji is set as representative of 
the first group. We choose thresholds thr 1 < £hr2, for 
example thrl = 20°, thr2 = 160°. For j m j 2 , • • • ,jn c we 
consider the angle a between tij 1 and tij. If a < thrl , 
the node j belongs to the first group. If a > thr2 and 
if the second group is empty, the node j becomes the 
representative node of the second group. If the second 
group is not empty, we consider the angle (3 between tij 
and where ko E • • • ,jn c } is the representative of 
the second group. If /3 < thrl , j is added to the second 
group. 

For the next search we replace tij 1 and tik 0 by the 
average normal vectors h\ and H 2 of group 1 and group 
2, respectively, and re-consider the nodes which could not 
have been assigned to one group in the first step. If the 
angle between tij and ni or tij and fi 2 is smaller than 
thrl, the node Xj is added to the corresponding group. 

If group 2 is empty, or if one of the groups consists 
of only a small number of nodes (e.g. < 5% of n c ), we 
start again by using another node as representative for 
the first group (e.g. j = j'2), since the node j = j\ 
could be an outlier. If no start node can be found, such 
that there exist two groups of nodes as described above, 
no topology change takes place. This can happen, if all 
weighted normals point in nearly the same direction. 

The method such provides that a topology change like 
splitting is not wrongly detected at locations where sev¬ 
eral nodes are just close to each other but their nor¬ 
mal vectors point in the same direction. From a mesh 
quality point of view, such meshes should be avoided; 
the nodes should be distributed equally over the surface. 
However, the detection of topology changes should be ro¬ 
bust enough not to wrongly detect a splitting at locations 
where there is just a high density of nodes. 

If both groups have a sufficient number of nodes, we 
proceed by considering the remaining normal vectors 
which could have not been assigned to one of the two 
groups. We again consider the angle between tij and 
Hi and tij and H 2 . If both angles are > thr3 (e.g. 
thr3 = 40°), the normal tij points in a complete dif¬ 
ferent direction compared to Hi and H 2 . Let No be the 
number of such points. If Nq exceed a predefined num¬ 
ber (e.g. 1/3 n c ), then there are more than two main 

groups of directions. In this case, there is a splitting or 
decrease of genus. Splitting and decrease of genus are 
handled similarly (see below for details how to modify 
the triangulations). Triangles close to the detected cube 
are deleted. If the remaining triangulation of the former 
surface consists of two connected components, a split¬ 


ting occurs. If the remaining triangulation is connected, 
a decrease of genus occurs. 

If No is zero or if it does not exceed the predefined 
number, possible remaining normal vectors are only sin¬ 
gle outliers and there are only two main groups of nor¬ 
mal vectors with nearly opposed normal vector. In this 
case, there is a merging or increase of genus. If there 
are nodes belonging to two different surfaces, a merging 
occurs. Otherwise, an increase of genus occurs. 

An illustration of the algorithm to detect and identify 
topology changes is given in Figure [l] 

3.3.3 Algorithm for Splitting and Decreasing of 
Genus 

We propose the following algorithm for a possible mod¬ 
ification of the surface triangulation after the detection 
of a splitting or decrease of genus. 

• Preparation and deletion of simplices: In case 
of splitting or decreasing genus, we consider the set 
of affected nodes X/, j E -S', which are located in 
the cube or in a neighbor cube, where the topology 
change has been detected. Let pe be the mean of the 
points {Xj : j E S'}. We delete all simplices with at 
least one vertex belonging to the set {Xj : j E S}. 
When deleting one simplex, we change the neigh¬ 
bor information of neighbor simplices at the corre¬ 
sponding edges to —1 (free edges). Deletion cre¬ 
ates two temporary holes in the surface(s). Sim¬ 
plices with two or three free edges are deleted as 
well, see Figure [2] As a result we either have two 
sets of connected simplices (—>> splitting) or one set 
of connected simplices (—» decrease of genus). 

• Set surface index (splitting only): The remain¬ 
ing simplices form two connected sets. For one set, 
we need to re-set the surface index. Let i be the 
surface index of the original surface. By splitting 
the total number of surfaces is increased to Ns + 1. 
Starting with one simplex with a free edge, we re-set 
its surface index from i to Ns + 1. Then, we con¬ 
sider its neighbor simplices and assign them to sur¬ 
face Ns + 1 also. By this procedure, the simplices of 
one of the two connected components are assigned 
one by one to the surface Ns + 1 by heritage, i.e. by 
use of neighbor information. 

• Generate new simplices: We close each of the 
two intermediate holes (where simplices have been 
deleted) by constructing new simplices at edges of 
simplices where the neighbor information has been 
set to —1. First we create two new points, each with 
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Figure 1: Illustration of the detection and identification of topology changes of surfaces. 
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Figure 2: Left: Part of the surface near a temporary hole. 
Free edges are drawn with red color. Right: Surface near 
the hole after deletion of simplices with more than one 
free edge. 




Figure 3: Improving mesh quality after a splitting or 
decrease of genus - Decrease the number of edges at the 
node pe- 


coordinates p#, one for each intermediate hole. In 
the next time steps the two nodes can move away 
from each other. If a is a simplex with a free edge 
given by X a j 1 and X a j 2 with no neighbor simplex, a 
new simplex is generated given by the vertices X a j 1 , 
X a j 2 and one of the two new nodes at pe- The new 
simplex inherits the surface index from a. 

• Improve mesh quality: By simply connecting all 
free edges with one of the two new vertices at pe, 
the two vertices can belong to a big number of sim- 
plices. Let a\ and a 2 be two of the newly generated 
simplices with one common edge. We construct 4 
new simplices from o\ and <72 such that the new 
vertex belongs to only one of the 4 new simplices, 
see Figure [3j Therefore the number of elements to 
which the vertices at pe belong to is approximately 
halved. We repeat this procedure until each of the 
newly created nodes at pe belongs to < 8 elements. 


3.3.4 Algorithm for Merging and Increasing of 
Genus 

We propose the following algorithm for a possible mod¬ 
ification of the surface triangulation after the detection 
of a merging or increase of genus. 


• Preparation and deletion of simplices: In case 
of merging or increasing genus, we consider the set 
of affected nodes Xj, j G S, which are located in 
the cube or in a neighbor cube, where the topol¬ 
ogy change has been detected. We delete all sim- 
plices with at least one vertex belonging to the set 
{Xj : j G S}. This generates temporary free edges, 
i.e. simplices exist which do not have a neighbor sim¬ 
plex at that edge. The neighbor index corresponding 
to the free edge is set to —1. This creates two in¬ 
termediate holes. Simplices with two or three free 
edges are deleted as well, see Figure [2] 

• Matching free nodes/edges: There exist now two 
sets of connected free edges. Let Ifree,k, k = 1, 2, be 
the set of nodes corresponding to the free edges (end 
points of the edges). We try to match the nodes of 
I free, l with the nodes of I free ,2 using the Hungarian 
method 2l| which is a combinatorial optimization 
algorithm. The Euclidean distance is used as cost 
criterion for matching two nodes. Since the number 
of nodes of the two sets need not be equal, there 
can be nodes which could not have been matched 
in the first step, see Figure [4] (top). Therefore, new 
nodes are created by bisection of simplices at a free 
edge. Finally, each node can be matched, see Fig¬ 
ure [4] (bottom). 

• Point/Edge merging: Since two matched nodes 
can have slightly different coordinates, they are re¬ 
placed by one node in the middle of the line con¬ 
necting the two nodes. The free edges of the two 
open holes are merged by identifying the matched 
nodes. The simplices, nodes and edges administra¬ 
tion needs to be adapted. The nodes in I free, 1 are 
updated; each point is replaced by the mid point 
between it and its matching partner. The nodes be¬ 
longing to I free, 2 are deleted. Half of the free edges 
are deleted. The edge, node and neighbor informa¬ 
tion of the simplices at the former holes need to be 
adapted. In case of merging, the surface index of all 
simplices belonging to the second surface is set to 
the surface index of the first surface. 


Note, that local refinement after a merging is typically 
necessary. This is automatically done, by a refinement 
method described in Section |3.4| If the surface grows 
locally near the former merging part, the simplices will 
become greater compared to the average simplex of the 
surface. In this case, the large simplices will be refined. 

The idea of creating two intermediate open holes and 
merging the two surfaces there is based on 12 . There 


however, each hole is restricted to consist of exactly four 
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Figure 4: A subset of the free nodes. Top: Intermedi¬ 
ate matching (red lines mark matching pairs). Bottom: 
Complete matching after inserting new nodes (red and 
green lines mark matching pairs). 

free edges. New triangles between the free edges are cre¬ 
ated instead of merging the edges. 

In our method, the seeking for close points (and there¬ 
fore close edges/simplices) is very efficient, since we make 
use of a background grid motivated by the method pre¬ 
sented in 26 . We extended this method originally in¬ 
tended for curves in the plane to topology changes of 
surfaces. We allow for intermediate holes with an ar¬ 
bitrary number of free edges. The hole size is of the 
magnitude of the grid size. 

3.4 Additional Computational Aspects 

3.4.1 Computations of regions and coefficients 

The computation of regions f }™ and coefficients c] 22 , the 
mean of uo in is done as follows: We assign each 
voxel of the three-dimensional image domain to a phase 
ft™. If a voxel is truncated by a surface, it is assigned 
to the phase to which the largest part belongs or to any 
of the two regions in case of two equal parts. Let S™ be 
the set of n™ voxels belonging to Q™. Then the approx¬ 
imation cf is set to 

rmi 

cr*= E «oU- ( 2 7) 

^ voxES™ 

The entire image domain needs to be considered only 
for mm 0. For m > 0, we only locally update the re¬ 
gions and re-compute the coefficients on this basis. For 
that, we consider a small band/tube of voxels around 
the current surfaces and look for changes of the region 
assignment. 

As the normal z/- 72 points from to the 

voxels close to the surface T™ can be assigned to the 


phase k + (i) or k~(i ), respectively. 

In the update step, we first set n™ = nj 72-1 and CJ! 1 = 
for k = 1,..., N r . For i = 1,..., Ns, all voxels in 
an environment of T- 72 are subsequently considered. Let 
a voxel vox be assigned to phase k G {fc + (i), k~(i)} and 
let l ^ k be the former phase index of the voxel. Then, 
we set 

< = < + 1, < = <-1, 

cr = cr+uoU, Q m = cr - u 0 \ vox . ( 28 ) 

After having considered all voxels close to the sur¬ 
faces, the coefficients are set to c™ = C™/n™ for k = 
1,... ,N r . 

3.4.2 Time Step Control 

We use a certain adaptive time step setting to control the 
speed of the evolution of the surface(s). Let At = r m de¬ 
note the (possibly variable) time step size. The time step 
size is controlled as follows: Let 5X™ in > 0, EA™ ax > 0 
with 5X™ in < EA™ ax be user-defined tolerances for the 
absolute value of the position difference in normal direc¬ 
tion. Let At > 0 be an initial time step size for m = 0 or 
the time step size of the previous time step for m > 0. 

We propose the following time step size control: 
Choose a factor X t G N (for example X t = 2 or X t = 10). 

1. Solve equation ( |26b| ) and set SX ™ +1 to the max¬ 
imum of | EA 772+1 . oj™j\ for i = 1,..., Ns and j 

1 . Niy. 

2. If 5X™ +1 > SX™ ax , set At to ±At and repeat 
step (i). 

3. Otherwise, if SX™ +1 < SX™ in , set At to X t At and 
repeat step (i). 

4. Otherwise, proceed by checking for topology changes 
(see above) and go to the next time step, i.e. set m 
to m + 1. 

The effect of this time step size control is simple: If there 
are too high changes in the position of the nodes in nor¬ 
mal direction (i.e. if the normal velocity is too high), the 
time step size will be decreased. This occurs if the sum 
of weighted curvature and external term is high. If the 
change in the position in normal direction is too small, 
the time step size will be increased to speed up the image 
segmentation process. 

3.4.3 Mesh Quality Aspects 

During the evolution of surface, it may be necessary to 
control the mesh quality. For example, if a surface con¬ 
tinuously grows, the simplices become larger and should 
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be refined if their area exceeds a certain threshold. Sim¬ 
ilarly, too small simplices should be deleted. 

For computing the matrix entries, cf. (23), we already 
need to compute the area of each simplex of the triangula¬ 
tion of the surface T™, i G {1,..., Ns }• Let ^desired > 0 
be a predefined desired area for one simplex. Let a > 0 
be a given factor (e.g. a = 2 or a = 10). If the area of 
a simplex exceeds a^desired? it will be refined by bisec¬ 
tion of its largest edge. Its neighbor simplex across the 
refinement edge will be also refined such that no hanging 
nodes remain. 

Local refinement is necessary to avoid too large sim¬ 
plices. Further, a mesh can also be continuously refined, 
for example when one starts with a small surface which 
globally grows. The triangles of the growing surface are 
refined one by one. Furthermore, the triangles which 
have one very large angle, i.e. an angle larger than a 
given threshold (e.g. > 160°), are also refined. 

If a simplex area is smaller than a certain percentage 
of the desired area AiesirecU for example smaller than 1%, 
the simplex is deleted. Further, it is also deleted if one of 
its three inner angles is smaller than a given threshold, for 
example smaller than 2°. When a simplex is marked for 
deletion, one or more neighbor simplices are also deleted. 

Mesh operations like deletion of triangles are rarely 
necessary. These operations are usually performed only 
a few times, for example close before or after topology 
changes. 

Figure [5] illustrates examples how the triangulation is 
adapted close to simplices which are marked for refine¬ 
ment or deletion. 

Our numerical method for surface evolution for image 
segmentation tasks is based on a numerical method de¬ 
veloped in [8] . This method provides a good mesh quality 
in many cases. However, if, for example, a pinch-off oc¬ 
curs, the mesh can get distorted and some routine for 
keeping a good mesh quality is needed. 

The idea of a mesh regularization method proposed 
in ( 9 ] is to induce or reduce the tangential motion of 
nodes along a surface. We use this method to control the 
tangential motion of nodes of surfaces during 3D image 
segmentation. For details, we refer to [9 . The system 
is replaced by a scheme which controls also the tan¬ 
gential motion of the nodes. 




Figure 5: Top: Refinement of a simplex (marked in gray). 
The neighbor simplex is also refined to avoid hanging 
nodes. Center: Deletion of a simplex with a too small 
area (marked in gray). Three neighbor simplices are also 
deleted. Bottom: Deletion of a simplex with a too small 
angle (marked in gray). One neighbor simplex is also 
deleted. 


i = 1 ,... , Ns , j = 1 , ..., W,v, perform the following 
steps for m = 0,1,..., M — 1: 


1. Compute the regions and the coefficients c™ 


3. 


k = 1,..., Nr, as described in Section 3.4.1 


2. Compute b m as defined in ( |24[ ) by using the coeffi¬ 
cients c™ of step[l] Compute X m+1 = + 4X m+1 

by solving the linear equation (26b). 


Check whether the time step size needs to be in¬ 
creased or decreased, see Section |3.4.2| If the time 
step size needs to be changed, repeat step [2] with the 
new time step size. 


4. Check whether topology changes occur and execute 
the topology change, see Section |T3j 

5. If necessary, refine too large simplices or delete too 
small simplices of the triangulation as described in 
Section 13.4.31 


3.5 Summary of the Image Segmentation 
Algorithm 

In summing up, we propose the following algorithm 
for segmentation of 3D images. Given a set of trian¬ 
gulated surfaces T 0 = (T?,..., T^- s ) and nodes X^ -, 


4 Results 

4.1 Artificial Test Images 

In this section we demonstrate the developed method for 
segmentation of 3D images with parametric active sur- 
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Figure 6: Splitting of a surface during 3D image segmen¬ 
tation. Surface(s) at step m = 0,100, 215, 300 at time 
t m = 0,0.01,0.0215,0.0298. 



Figure 7: Merging demonstration. Surface(s) at step 
m = 0, 50,80, 200 at time t m = 0, 0.005, 0.008, 0.02. 


faces. We first study four examples of artificial test im¬ 
ages to demonstrate the ability of the method to detect 
different topology changes (splitting, merging, increase 
and decrease of genus). 

In the first experiment, we demonstrate how a surface 
is split in two surfaces. We consider an artificial image 
defined on an image domain given given by the cuboid 
n = [-2.5, 2.5] x [-1.5,1.5] x [-1.5,1.5]. The image in¬ 
tensity function is defined by 


of two surfaces to one single surface. The initial surfaces 
in the next example are two balls. The image domain is 
given by = [—1.2,1.2] x [—0.8, 0.8] x [—0.8, 0.8] and the 
image intensity function is defined by 

f 0 if ||x|| < 0.6, 
uq : Q — )> M, uq(x) = < 

I 1 else. 


uq : Q —>> M, uq(x) = < 


0 if \\x- (-1.2,0,0) T || < 0. 
V||f- (1.2,0,0) T || < 0.8, 

1 else. 


The three-dimensional image contains two balls centered 
at (±1.2,0,0 ) t G M 3 with radius 0.8. The segmentation 
process is started using a cylinder-like surface as initial 
surface placed in the center of the cuboid. Figure [6] shows 
the surface at different time steps. 

For weighting the curvature term and the forcing term 
for the image segmentation, the parameters a = 1 and 
A = 100 are used. At time step m = 214 a splitting 
of the evolving surface occurs. In the subsequent itera¬ 
tions steps, the two new surfaces each evolve to a ball. 
To detect the topology change, we use an auxiliary back¬ 
ground grid with grid size a = 0.025 as described in Sec¬ 
tion |3.3.1| A cube of the grid is considered for possible 
topology changes if more than A^etect = 10 nodes are 
located inside the cube. We further use the parameters 
thrl = 30°, thr2 = 150° and thr3 = 40°. 

Furth ermor e, we also perform a time step control (cf. 
Section 3.4.2) using the thresholds 5X™ in = 0.003 and 


5X™ ax = 0.05. For almost each time step, a time step 
size of At = 10 -4 was used. Only immediately after the 
splitting, for two time steps, At was reduced to 10 -5 to 
avoid a too fast retraction of the newly generated surfaces 
close to the former splitting point. 

The reversed topology change of splitting is a merging 


As weighting parameters a = 2 and A = 60 are used. 
Time step control is performed applying the thresholds 
5X™ in = 0.001 and SX™ ax = 0.02. No change of the 
time step size is necessary in this example; the time step 
size At = 10 -4 need not be changed throughout the evo¬ 
lution. To detect the merging, a = 0.03, Ahetect = 10 
and thrl = 20°, thr2 = 150° and thr3 = 40° are used. 
The resulting surfaces of this experiment at different time 
steps are shown in Figure [7| 

Since the surface grows continuously, some simplices 
have to be refined as described in Section 13.4.31 The 
desired area for one simplex is Adesired = 0.001; a simplex 
is refined by bisection of its largest angle if its area is 
larger than a certain factor of Ad es ired- A simplex is also 
bisected if one angle is larger than 170°. A simplex is 
deleted if one angle is smaller than 2° or if its area is 
smaller than 1% of the desired area of Adesired- 

In the next examples, we demonstrate another kind of 
topology changes: increase and decrease of the genus of 
a surface. Therefore, we consider an image segmentation 
example where a sphere should evolve to a torus. The 
image intensity function is given by 



if (y/x\ + x\ - R) 2 + x\ < r 2 , 
else, 


(29) 


where R = 1.2 and r = 0.4 are used here. 

Figure [8] shows the surface, its mesh (cross-section of 
the mesh) at different time steps. For this example we 
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Figure 8: Demonstration of an increase of the genus of 
a surface. Surface, mesh and cross-section at step m = 
0,325,500 (row-wise) at time t m = 0,0.325,0.5. Column 
1-2: surface and mesh (cross-section). 
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Figure 10: Time step sizes during the evolution of the 
torus to a ball. 


apply a = 1, A = 60 (weighting parameters). The topol¬ 
ogy change is detected using a = 0.0565, A detec t = 8 and 
thrl = 20°, thr2 = 150° and thr3 = 40°. As parameters 
to control the refinement, the desired triangle area is set 
to A desired = 0.005, and the angles 170° and 2° are used 
for bisection or deletion of a triangle, respectively. 

Further, 5X™ in = 0.01 and 5X™ ax = 0.1 are applied 
as thresholds for the time step size control. Throughout 
the evolution, there was no need to change the initial 
time-step size of At = 10 -3 . 

At time step m = 325, two different parts (top and 
bottom) of the surface with nearly opposite normal vec¬ 
tor nearly touch. A topology change is detected and a 
small hole occurs. The genus of the surface is increased 
from g = 0 to g = 1. At time step m = 500, the 3D 
object, a torus, is detected; its boundary is represented 
by the surface. 

Finally, we present an example where a torus is used 
as initial surface and a sphere should be detected. The 
image intensity function is given by 


Uq : 0, —» M, u 0 (x) 


0 if \\x || < 0.8, 
1 else. 


% • 


Figure 9: 3D image segmentation example where a torus 
evolves to a ball. Surface at step m = 0,425, 500,1000 
(row-wise) at time t m = 0,0.0425,0.04406,0.27316. Col¬ 
umn 1-2: surface and mesh (cross-section). 


Figure [9] shows the surface at several time steps. As 
weighting parameters a = 1 and A = 20 are applied. For 
the detection of the decrease of genus, the parameters 
a = 0.025, Adetect = 20 and thrl = 20°, thr2 = 150° and 
thr3 = 40° are used. 


The time step size is controlled using the thresholds 
5X™ in = 0.0005 and SX™ ax = 0.01. Figure [To] shows the 
time step sizes during the image segmentation process. 
After the topology change the time step size is decreased 
from 10 -4 to 10 -5 . Later it is increased to speed up the 
segmentation. 
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Figure 11: Lung segmentation: Surfaces (row 1) and 
cross-sections (row 2: z = 80,150,200, row 3: y = 
80,150, 200) at m = 0 at time t = 0. The original images 
are from the Lung Image Database Consortium image 
collection (LIDC-IDRI) of The Cancer Imaging Archive 
(TCIA), see [31], |], [30 .. 


4.2 Segmentation of Medical 3D Images 

In this section, we apply the segmentation method for 
three-dimensional images to medical image data. Seg¬ 
mentation of medical images is a challenging task due to 


possible high noise and image artifacts, see 35 


3D image data often consists of a set of 2D slice im¬ 
ages generated by radiology scans, for example computed 
tomography (CT) and magnetic resonance (MR) scans. 
With a 3D image segmentation technique, one can seg¬ 
ment organs (heart, lung, abdomen, liver, etc.) or tumors 
from their environment. The output, i.e. the resulting 
surface, serves as a reconstruction and visualization of 
the medical object and could be used for further medical 
analysis and diagnostic purposes: After the segmenta¬ 
tion, one can compute the area of the triangulated sur¬ 
faces and the volume of the enclosed regions. The area 
of the surfaces and the volume of the regions could be 
used for example to analyze if a tumor has been growing 
in the time between two radiological examinations. 

First, we consider a sample 3D image of the 
Lung Image Database Consortium image collec¬ 
tion (LIDC-IDRI) of The Cancer Imaging Archive 
(TCIA) (https://wiki.cancerimagingarchive.net/ 


display/Public/LIDC-IDRI), see 31 


mgarchive. 

> iiifl 


The 


1 The author acknowledges the National Cancer Institute and 
the Foundation for the National Institutes of Health, and their 




Figure 12: Lung segmentation: Surfaces (row 1) and 
cross-sections (row 2: z = 80,150,200, row 3: y = 
80,150, 200) at m = 100 at time t = 10. The original 
images are from the Lung Image Database Consortium 
image collection (LIDC-IDRI) of The Cancer Imaging 
Archive (TCIA), see [3^, [|, [30]. 




Figure 13: Lung segmentation: Surfaces (row 1) and 
cross-sections (row 2: z = 80,150,200, row 3: y = 
80,150, 200) at m — 600 at time t = 60. The original 
images are from the Lung Image Database Consortium 
image collection (LIDC-IDRI) of The Cancer Imaging 
Archive (TCIA), see (31 , |], [30]. 
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data set consists of diagnostic CT scans. The original 
data set consists of 2D slice images stored as DICOM 
files. The files are first preprocessed to cuboid 3D images 
with N x x N y x N z voxels, here: N x = 445, N y = 310 


and N z = 250. 

Figures PPi show the evolving 3D surfaces and six 
representative 2D cross-sections at different time steps 
m = 0,100,600. In the subfigures showing 2D cross- 
sections, the image cross-sections for constant z (in de¬ 
tail z = 80,150,200) and constant y (in detail y = 
80,150, 200) are drawn as well as the intersection points 
of the surfaces’ edges with the cross-section planes. 

For the image segmentation the weight of the curvature 
term is set to a = 10, the weight of the external forcing 
term to A = 1000. As parameters for the time step size 
control, 5X™ m = 0.05 and SX™ ax = 2 are used. The 
time step size At = 0.1 need not be changed during the 
segmentation. 

In a second experiment, we consider an experiment 
where a topology change occurs. We perform again a 
lung segmentation starting now with one initial surface 
which is split into two surfaces. Figure [l4p8] show the 
surface(s) at time step m — 0, 50,100, 500,900 as well 
as cross-sections of the image and of the surface(s). For 
the cross-sections, we consider the planes given by z = 
80,120,160 and y = 50, 64, 80. 

The splitting occurs at time step m — 50. To detect 
the topology change, we use an auxiliary background grid 
with grid size a = 2. A cube of the grid is considered for 
possible topology changes if more than detect = 8 nodes 
are located inside the cube. Further, we use the param¬ 
eters thrl = 30°, thr2 = 150° and thr3 = 40°, recall 
Section |~3.3.2| After the splitting, the two surfaces grow 
and new triangles are created by bisection of too large 
triangles. For the segmentation we use the parameters 
(7 = 1 and A = 20. The time step size is set to At = 0.2 
with time step control using = 2, 5X™ in = 0.1. 

However, no increase or decrease of the time step size is 
necessary. 

As postprocessing step, we compute the volume of the 
two enclosed regions and the area of the region bound¬ 
aries. The right lung of the patient, i.e. the left surface 
in the Figure 18, has an area of A\ = 3.309 • 10 4 and a 
volume of V\ = 2.691 • 10 5 (CT images are mirror im¬ 
ages). The left lung of the patient (right surface in the 
figure) has an area of A 2 = 2.801 • 10 4 and a volume of 
V 2 = 1.923 • 10 5 . Thus, as expected, the volume of the 
right lung is larger compared to the left lung. Note, that 
we handle a voxel as a cube with side length 1, resulting 
in values of magnitude 10 4 for the area and 10 5 for the 


critical role in the creation of the free publicly available LIDC/IDRI 
Database. 



20 70 120 20 70 120 20 70 120 


Figure 14: Lung segmentation with splitting: Surface at 
different viewing angles (row 1) and cross-sections (row 2: 
z = 80,120,160, row 3: y = 50,64, 80) at m = 0 at time 
t = 0. Credits (original CT images): C. Stroszczynski, 
Radiology, University Hospital Regensburg. 
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Figure 15: Lung segmentation with splitting: Surfaces 
(row 1) and cross-sections (row 2: z = 80,120,160, row 
3: y = 50, 64, 80) at m = 50 at time t = 10. Credits 
(original CT images): C. Stroszczynski, Radiology, Uni¬ 
versity Hospital Regensburg. 


Figure 16: Lung segmentation with splitting: Surfaces 
(row 1) and cross-sections (row 2: z = 80,120,160, row 
3: y = 50,64, 80) at m = 100 at time t = 20. Credits 
(original CT images): C. Stroszczynski, Radiology, Uni¬ 
versity Hospital Regensburg. 
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Figure 17: Lung segmentation with splitting: Surfaces 
(row 1) and cross-sections (row 2: z = 80,120,160, row 
3: y = 50, 64, 80) at m = 500 at time t = 100. Cred¬ 
its (original CT images): C. Stroszczynski, Radiology, 
University Hospital Regensburg. 


Figure 18: Lung segmentation with splitting: Surfaces 
(row 1) and cross-sections (row 2: z = 80,120,160, row 
3: y = 50,64, 80) at m = 900 at time t = 180. Cred¬ 
its (original CT images): C. Stroszczynski, Radiology, 
University Hospital Regensburg. 
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volume. If details on the acquisition system of the CT 
images are known (like the slice thickness, and the height 
and width of one pixel of a slice image), the area and the 
volume can be computed precisely and can be expressed 
in the metric system for practical interpretation of the 
values. 

5 Conclusion 

We presented a new parametric method for segmentation 
of 3D images. We considered extensions of the Mumford- 
Shah and Chan-Vese functional for 3D image segmenta¬ 
tion by active surface. For the time-dependent surfaces, 
we proposed a parametric scheme and introduced an effi¬ 
cient numerical scheme based on a finite element approx¬ 
imation. A novel method to detect and perform topology 
changes of the surfaces has been presented which uses a 
virtual auxiliary background grid. Due to the fact that 
for the main computations only a two-dimensional grid is 
used, the developed method is very efficient from a com¬ 
putational point of view. Several artificial images have 
been studied to demonstrate splitting and merging of sur¬ 
faces, and increase and decrease of the genus of a surface. 
We successfully applied our method to real medical 3D 
image data from computed tomography, including an ex¬ 
ample with a topology change. 
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